#cd /Paterson/Datafiles/grouse/mapping
#grep -n '>' grouse1.tsv |sed 's/:>/ /' > grouse1_lines.txt
#grep -n '>' grouse2.tsv |sed 's/:>/ /' > grouse2_lines.txt
#grep -n '>' grouse3.tsv |sed 's/:>/ /' > grouse3_lines.txt


readTsvFile <- function(tsv.file,line.file){
	#Function to get reference points for contigs within tsv file
	#tsv.file is the name of the tsv file
	#line.file is created by one of the grep commands above
	lines.df <- read.table(line.file,col.names=c("line","contig","start"),stringsAsFactors=FALSE)
	lines.df$start <- lines.df$line+1
	lines.df$stop <- 0
	lines.df$stop[-length(lines.df$stop)] <- 
	lines.df$line[-1] -1
	tmp <- system(paste('wc -l',tsv.file),intern=TRUE)
	lines.df$stop[length(lines.df$stop)] <- strsplit(tmp," ")[[1]][2]
	lines.df
	}

getContigFromTsv <- function(contig,lines.df,tsv.file){
	#extract info for a contig from tsv file 
	#lines.df created by readTsvFile
	
		
	tsv.df <- data.frame(position=numeric(0),reference=character(0),consensus=character(0),quality.score=integer(0),depth=integer(0),signal=numeric(0),std.deviation=numeric(0))
	
	
	tmp <- lines.df[lines.df$contig==contig,]
	if(nrow(tmp)==0) return(tsv.df)
	
	for(i in 1:nrow(tmp)){
		tmp.cmd <- paste("sed -n \'",tmp$start[i],",",tmp$stop[i],"p\' ",tsv.file,"> tmp_contig.txt",sep="")
		system(tmp.cmd)
		tmp.contig <- read.table('tmp_contig.txt',stringsAsFactors=FALSE)
		tsv.df <- rbind(tsv.df,tmp.contig)
		}
	names(tsv.df) <- c("position","reference","consensus","quality.score","depth","signal","std.deviation")
	tsv.df
	}


grouse1_lines <- readTsvFile('mapping/grouse1.tsv','mapping/grouse1_lines.txt')
grouse2_lines <- readTsvFile('mapping/grouse2.tsv','mapping/grouse2_lines.txt')
grouse3_lines <- readTsvFile('mapping/grouse3.tsv','mapping/grouse3_lines.txt')


#write functions to produce pretty plots of contigs, coverage and snps
embl.file <- "/Paterson/Datafiles/grouse/embl_exon2/cDNA_183-1_exon.embl"

plot.contig.coverage <- function(embl.file,snp.df,grouse_lines,tsv.files,...){
	#plots coverage, snps and cds for a sequence captured gene
	#grouse lines is a list of data.frames produced by readTsvFile
	#tsv.files is a character vector giving names of tsv files

	embl.lines <- readLines(embl.file)
	embl.CDS.lines <- grep('^FT[ ]*CDS',embl.lines)
	
	embl.product.lines <- grep('/product',embl.lines)
	embl.CDS <- rep("x",length(embl.product.lines))
	for(cdsi in 1:length(embl.product.lines)){
		tmp <- embl.lines[embl.CDS.lines[cdsi]:(embl.product.lines[cdsi]-1)]
		tmp <- sub('^FT[ ]*','',tmp)
		tmp <- sub('^CDS[ ]*','',tmp)
		embl.CDS[cdsi] <- paste(tmp,collapse="")
		}
	#embl.CDS <- gsub('FT[ ]*','',embl.CDS)
	#embl.CDS <- sub('CDS[ ]*','',embl.CDS)
	embl.CDS.dir <- "forward"
	if(substr(embl.CDS,1,4)=="comp") embl.CDS.dir <- "reverse"
	embl.CDS <-sub('[a-z]*\\(','',embl.CDS)
	embl.CDS <- sub('\\)$','',embl.CDS)

	#create a complex list
	embl.CDS.list <- lapply(strsplit(embl.CDS,','),strsplit,'\\.\\.')
	embl.CDS.start <- lapply(embl.CDS.list,function(X){sapply(X,function(X){as.numeric(X[1])})})
	embl.CDS.stop <- lapply(embl.CDS.list,function(X){sapply(X,function(X){as.numeric(X[2])})})

	embl.length <- as.numeric(sub('^FT[ ]*source[ ]*1\\.\\.','',embl.lines[grep('^FT[ ]*source',embl.lines)]))

	
	contig.names <- strsplit(paste(sub('^CC[ ]*','',embl.lines[grep('^CC',embl.lines)]),collapse=" "),'; ')[[1]]

	#where are the contigs
	contig.lines <- grep('/note=\"contig',embl.lines)
	contig.direction <- 	factor(rep("forward",length(contig.lines)),levels=c("forward","reverse"))
	contig.dir.txt <- embl.lines[contig.lines+1]
	contig.direction[grep('reverse',contig.dir.txt)] <- "reverse"
	
	contig.pos <- sub('^FT[ ]*misc_feature[ ]*','',embl.lines[contig.lines-1])
	contig.mat <- sapply(contig.pos,function(X){as.numeric(strsplit(X,'\\.\\.')[[1]])})
	#1st row has start points, 2nd row has end points


	#get blast alignment, ie what went onto the array
	blast.lines <- grep('/note=\"blast hit to query',embl.lines)
	blast.pos <- sub('^FT[ ]*misc_feature[ ]*','',embl.lines[blast.lines-1])
	blast.mat <- sapply(blast.pos,function(X){as.numeric(strsplit(X,'\\.\\.')[[1]])})

	#get the depth of coverage ready for plotting
	contig.tsv <- vector('list',length=3)
	names(contig.tsv) <- c('grouse1','grouse2','grouse3')
	max.depth <- 1

	for(cdsi in 1:length(contig.names)){
		contig.tsv[[1]][[cdsi]] <- getContigFromTsv(contig.names[cdsi],grouse_lines[[1]],tsv.files[1])
		contig.tsv[[2]][[cdsi]] <- getContigFromTsv(contig.names[cdsi],grouse_lines[[2]],tsv.files[2])
		contig.tsv[[3]][[cdsi]] <- getContigFromTsv(contig.names[cdsi],grouse_lines[[3]],tsv.files[3])
		depth.cdsi <- max(c(contig.tsv[[1]][[cdsi]]$depth,contig.tsv[[2]][[cdsi]]$depth,contig.tsv[[3]][[cdsi]]$depth))
		if(depth.cdsi > max.depth) max.depth <- depth.cdsi
	
	
		if(contig.direction[cdsi]=="forward"){
			contig.tsv[[1]][[cdsi]]$rel.pos <- contig.mat[1,cdsi] -1 + contig.tsv[[1]][[cdsi]]$position
			contig.tsv[[2]][[cdsi]]$rel.pos <- contig.mat[1,cdsi] -1 + contig.tsv[[2]][[cdsi]]$position
			contig.tsv[[3]][[cdsi]]$rel.pos <- contig.mat[1,cdsi] -1 + contig.tsv[[3]][[cdsi]]$position
			}else{ #reverse
			contig.tsv[[1]][[cdsi]]$rel.pos <- contig.mat[2,cdsi] - contig.tsv[[1]][[cdsi]]$position +1
			contig.tsv[[2]][[cdsi]]$rel.pos <- contig.mat[2,cdsi] - contig.tsv[[2]][[cdsi]]$position +1
			contig.tsv[[3]][[cdsi]]$rel.pos <- contig.mat[2,cdsi] - contig.tsv[[3]][[cdsi]]$position +1
			}
	
		}
	#quality control on SNPs
	snp.df2 <- snp.df[apply(snp.df[,grep('cov',names(snp.df))],1,function(X){sum(as.numeric(X))>30}),]

	#do the plots
	plot(x=0,y=0,type="n",xlim=c(0,embl.length*1.1),ylim=c(1.5,10),ylab="",xlab="base pairs")
	lines(x=c(1,embl.length),y=c(5,5))
	rect(xleft=contig.mat[1,],xright=contig.mat[2,],ybottom=5,ytop=5.5)
	rect(xleft=blast.mat[1,],xright=blast.mat[2,],ybottom=5,ytop=5.5,density=20)
	
	for(cdsi in 1:length(embl.CDS.start)){
		rect(xleft=embl.CDS.start[[cdsi]],xright=embl.CDS.stop[[cdsi]], ybottom=5+cdsi,ytop=5.5+cdsi)
		for(cdsj in 1:(length(embl.CDS.start[[cdsi]])-1)){
			lines(x=c(embl.CDS.stop[[cdsi]][cdsj],rep(mean(c(embl.CDS.stop[[cdsi]][cdsj],embl.CDS.start[[cdsi]][cdsj+1])),2),embl.CDS.start[[cdsi]][cdsj+1]),
			y=c(5.25+cdsi,5.25+cdsi+0.2,5.25+cdsi+0.2,5.25+cdsi))
			
			}
		}

	cov.pretty <- pretty(c(0,max.depth),n=3) 
	for(i in 1:length(contig.names)){
		points(contig.tsv[[1]][[i]]$rel.pos,2*contig.tsv[[1]][[i]]$depth/cov.pretty[length(cov.pretty)]+2,pch=".")
		points(contig.tsv[[2]][[i]]$rel.pos,2*contig.tsv[[2]][[i]]$depth/cov.pretty[length(cov.pretty)]+2,pch=".")
		points(contig.tsv[[3]][[i]]$rel.pos,2*contig.tsv[[3]][[i]]$depth/cov.pretty[length(cov.pretty)]+2,pch=".")

		}

	lines(x=rep(embl.length*1.05,2),y=c(2,4))
	lines(x=embl.length*c(1.05,1.06),y=c(2,2))
	lines(x=embl.length*c(1.05,1.06),y=c(3,3))
	lines(x=embl.length*c(1.05,1.06),y=c(4,4))
	text(x=embl.length*1.06,y=2,pos=4,labels="0")
	text(x=embl.length*1.06,y=3,pos=4,labels=max(cov.pretty)/2)
	text(x=embl.length*1.06,y=4,pos=4,labels=max(cov.pretty))


	points(x=snp.df2$start.pos[snp.df2$dn==FALSE], y=rep(4.5,length(snp.df2$start.pos[snp.df2$dn==FALSE])),pch="|")
	points(x=snp.df2$start.pos[snp.df2$dn==TRUE], y=rep(4,length(snp.df2$start.pos[snp.df2$dn==TRUE])),pch="|")


	}

#plot.contig.coverage("/Paterson/Datafiles/grouse/embl_exon2/cDNA_183-1_exon.embl",hcCoveragelist2[['cDNA_183-1']],grouse_lines=list(grouse1=grouse1_lines,grouse2=grouse2_lines,grouse3=grouse3_lines),tsv.files=c('mapping/grouse1.tsv','mapping/grouse2.tsv','mapping/grouse3.tsv'))

#plot.contig.coverage("/Paterson/Datafiles/grouse/embl_exon2/cDNA_469-1_exon.embl",hcCoveragelist2[['cDNA_469-1']],grouse_lines=list(grouse1=grouse1_lines,grouse2=grouse2_lines,grouse3=grouse3_lines),tsv.files=c('mapping/grouse1.tsv','mapping/grouse2.tsv','mapping/grouse3.tsv'))

#cDNA_1686-1 looks ok, though have to explain 5' UTR
#cDNA_1686-1	ENSGALP00000018553	TNF receptor-associated factor 3 (CD40 receptor-associated factor 1)(CRAF1)(CD40-binding protein)(CD40BP)(LMP1-associated protein)(LAP1)(CAP-1) [Source:UniProtKB/Swiss-Prot;Acc:Q13114]	sp|Q13114|TRAF3_HUMAN	Q13114	TNF receptor-associated factor 3



svg("grouse_analysis/cDNA_1686-1.svg")
plot.contig.coverage("/Paterson/Datafiles/grouse/embl_exon2/cDNA_1686-1_exon.embl",hcCoveragelist2[['cDNA_1686-1']],grouse_lines=list(grouse1=grouse1_lines,grouse2=grouse2_lines,grouse3=grouse3_lines),tsv.files=c('mapping/grouse1.tsv','mapping/grouse2.tsv','mapping/grouse3.tsv'))
dev.off()





